library(vegan)
speabu<-read.csv('abund.csv',header=TRUE,row.names=1)
speabu.log<-log(speabu[,-1]+1)
#perform NMBD on log transformed data distrance is distance coefficient to use (default is bray-curtis), k is the number of ordnation axes to generate, autotransform default is TRUE, will transform the data first (usually square root), trymax is the number of random starts to try, default is 20
spe.nmds<-metaMDS(speabu.log,distance='bray', k=2, autotransform=FALSE, trymax=100)
#look at it
spe.nmds
#get a list of objects resulting from the analysis
names(spe.nmds)
#stress is kind of high we could transform to presence absence or we could add a 3RD axis
spe.nmds2<-metaMDS(speabu.log,distance='bray', k=3, autotransform=FALSE, trymax=100)
#look at it
spe.nmds2
#evaluate the number of dimensions with a scree plot
nmds.scree(speabu.log, distance='bray', k=10, autotransform=FALSE, trymax=20)
#or try the randomization thing
nmds.monte(speabu.log, distance='bray', k=3, autotransform=FALSE, trymax=20)
#how good a job does NMDS do? check by looking at corrlation between the caclulated dissimilarities and the plotted values. Specifically, plot original dissimilarities and euclidean distances in the ordination using stressplot()
stressplot(spe.nmds2)
#look at a plot of the NMDS
plot(spe.nmds, type='n')
text(spe.nmds, labels=row.names(speabu))
#if we'd like to see how a particular descriptor (in this case, a weird fish abundance) changes with location we can make te symbol size proportional to log abundance
plot(spe.nmds, type='n')
points(spe.nmds,cex=speabu.log$ROSYDACE)
#if we want to look at the sample scores, which are the coordinates of the samples in k-dimensional ordination, they are stored in the object names points
spe.nmds$points
#to calculate species loadings (i.e. variable weights) on each derived axis, use envfit with the NMDS scores
vec.sp<-envfit(spe.nmds$points, speabu.log, perm=1000)
vec.sp
#NOW these loadings can be plotted
ordiplot(spe.nmds, choices=c(1,2), type='text', display='sites', xlab='Axis 1', ylab='Axis 2',)
plot(vec.sp, p.max=0.01, col='blue')